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ABSTRACT 

We present our simulation tool JCMmode for calculating propagating modes of an optical waveguide. As ansatz 
functions we use higher order, vectorial elements (Nedelec elements, edge elements). Further we construct 
transparent boundary conditions to deal with leaky modes even for problems with inhomogeneous exterior 
domains as for integrated hollow core Arrow waveguides. We have implemented an error estimator which steers 
the adaptive mesh refinement. This allows the precise computation of singularities near the metal's corner of a 
Plasmon-Polariton waveguide even for irregular shaped metal films on a standard personal computer. 

• i-H ' 

Keywords: Leaky Modes, Nano-Optics, Plasmon-Polariton Modes, Arrow Waveguide, Finite-Element-Method, 
Q ■ Pole Condition, PML 

1. INTRODUCTION 

The computation of propagating modes of an optical waveguide is one of the central tasks in the optical component 
i-C ' design. In mathematical modeling this corresponds to a quadratic eigenvalue problem in the sought propagation 
constant k z } Beyond "true" eigenmodes with finite energy in the cross section there exist so-called "leaky 
modes" which are solutions to Maxwell's equations but with typically increasing field intensity for a growing 
distance to the waveguide core. 2 ' 3 These leaky modes must satisfy a further asymptotic boundary condition 
for large distances to the waveguide core. Analog to scattering problems, one demands that there is no energy 
transport from infinity within the cross section. 4 ' 5 To bring this into a mathematical form, we split the 
cross section R 2 into a bounded interior domain Oi„t and an exterior domain f2 e xt, that is R 2 = J7i n t U fi C xt- 
. For a homogeneous exterior domain f2 ex t (with constant permittivity and permeability) the correct asymptotic 
boundary condition is the well known Silver-Miiller condition. 6 Uranus and Hoekstra use a BGT-like transparent 
boundary condition based on this asymptotic boundary condition. 3 Besides a poor convergence with the size of 
the computational domain, this asymptotic boundary condition is wrong for inhomogeneous exterior domains. 4 
But, many waveguide structures are composed of layers with an immense lateral expansion compared to the 
. — i ' waveguide core diameter. These structures are best modeled in the way that the layers reach infinity. To 
deal with such inhomogeneous exterior domains in a rigorous manner, Schmidt has proposed the pole condition 
concept for the definition of asymptotic boundary conditions. ' 5 We briefly introduce this concept in Section 3. 
O r Further we show the connection of this concept to a modified PML method proposed by the authors. 7 In the 
Section 5 we explain how to discretize the modified PML method and how to couple the transparent boundary 
condition with the interior finite element discretization. In the last section we demonstrate the ability of our 
method for challenging problems in modern optical waveguide design. 

Alternatively to the modified PML method Schmidt has presented a numerical approach which is directly 
based on the pole condition. The authors will compare these two methods in a succeeding paper. 
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2. LIGHT PROPAGATION IN A WAVEGUIDE 

Starting from Maxwell's equations in a medium without sources and free currents and assuming time-harmonic 
dependence with angular frequency u> > the electric and magnetic fields 

E(x, y, z, t) = E(x, y, z)e^\ U(x, y, z, t) = H(x, y, z)e~ i ^ t , 

must satisfy 

V x E = iufiH, V • eE = 0, 

V x H = -iweE, V • /iH = 0. 



Here e denotes the permittivity tensor and \i denotes the permeability tensor of the materials. In the following we 
drop the wiggles, so that E — > E, H — > H. From the equations above we then may derive (by direct substitution) 
the second order equation for the electric field 

Vx^VxE-w 2 fE = 0, 
V • eE = 0. 

A similar equation holds true for the magnetic field - one only must replace E by H and interchange e and [i. 
Observe that any solution to the first equation also meets the divergence condition (second equation). This is 
because V • Vx =0. 

To recover the underlying structure we rewrite these equations in differential form, 

di[i~ 1 die — uj 2 ee = 0, (la) 
d 2 ee = 0. (lb) 

A reader not familiar with this calcalus may replace the exterior derivatives do, d\, d\ with classical differential 
operators, d a — > V, d\ — > Vx and d 2 — > V-. Here, the electric field appears as a differential 1-form, e = 
e x dx + e y dy + e z dz, whereas the material tensors act - from a more mathematical point of view - as operators 

e, ii : Alt 1 -» Alt 2 . 

In order to derive a weak formulation we define the following function spaces on the domain ft = R 3 

Hl c = {0e Alt°|V0e (L 2 Xoc f} 

ffioc(curl) = {e G Alt 1 | (e x ,e y ,e z ) e {L 2 Xoc f , V x (e x , e y , e z ) T e (L? oc ) 3 } 
Fi oc (div) = {d e Alt 2 | (d x ,d y ,d z ) G (Lf oc f, V • (d x ,d y , d z f G Ll c ) 

The weak form to Equations (1) now reads 

(/i _1 c?ie A d\v — Lu 2 (ee) A v) = (2a) 



R 3 

(ee) A d Q p = (2b) 



for all v G iJi oc (curl) and p G H\ oc with compact support. 

An optical waveguide is an invariant structure in one spatial direction which we assume to be the z - direction 
of a cartesian coordinate system. A propagating mode is a solution to the above time-harmonic Maxwell's 
equations such that the electric field E depends harmonically on the spatial coordinate z, 



E(x,y,z)=E(x,y)e^- z . 
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Figure 1. Prismatoidal coordinate system. Each segment Qj is the image of a reference element under a bilinear mapping 
B l j OC . These local mappings are combined to a global mapping B which is continuous in n. 

Hence a propagating mode travels along the z-dircction. The scalar quantity k z is called propagation con- 
stant. Let us denote by H\ OCt k z (curl) the subspace of fields in i?i oc (curl) which depends on z as e(x,y, z) = 
e(x,y, z) exp(ik z z). The spaces H^ oc k and i?ioc,fcz(div) are defined accordingly, ft is sufficient to restrict the 
variational problem (2) on the cross section z — 0. As mentioned in the introduction a propagating mode should 
not only solve Maxwell's equations but should also transport no energy within the cross section from infinity 
that is it should be purely outgoing in the cross section. The precise definition of what purely outgoing means is 
given in the next section. The weak waveguide problem is summarized in the following Problem 1. 

Problem 1 (Weak Waveguide Problem). Find k z such that there exists a field e £ H\ oc ^ kz (curl) which 
is purely outgoing in the cross section and which satisfies 

(3a) 
(3b) 
for any v £ i?ioc,fc 2 (curl), p £ H^ c k ^ with compact support in x and y. 

3. LEAKY MODES AND OUTGOING BOUNDARY CONDITION 

We now address the definition purely outgoing in Problem 1. From a physical point of view, any propagating 
mode is admissible as long as there is no energy transport in the cross section from infinity. As mentioned in 
the introduction to this paper we want to define the transparent boundary condition with the help of the pole 
condition concept, 4 which we now detail for the one dimensional case. 

Let us assume that the permittivity and permeability are only dependent on x, e = e(x), p = /j,(x) and are 
constant in the right exterior domain I + = [0, +oo). Then a TE mode satisfies the Helmholtz equation 

-d xx E y {x) + k 2 z E y (x) - Lo 2 iicE y {x) = 0, x e /+ 

with general solution 



If we define the square root so that 3?y / u; 2 ^.e — k z > the first part is an outgoing wave and the second part is 
an incoming wave. Therefore, as "physical" boundary condition we must enforce B = 0. 



/ (fi 1 die A div — uj 2 (ee) A v) 
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Figure 2. Discretization of the interior domain and rays in the exterior domain (left picture). Geometry with represen- 
tation of the refractive index distribution (right picture). Infinite waveguide: &2 = 1.32, background: fci = 0.29. 



Let us regard the Laplace transform of E y , 

/>OC 

CE y = / E y (x)e~ sx dx = 
Jo 



A 



B 



s — i\J uS 1 [it — k 2 s + iy/ uj 2 iie — k 2 



We see that the incoming wave produces a pole at s = —i^u; 2 ^ — k 2 . Hence B = is equivalent to the fact 
that the Laplace transform of the solution is holomorphic in the lower complex half plane. This is precisely the 
pole condition for the one dimensional case: 

A solution to Helmholtz equation (4) is purely outgoing if its Laplace transform is holomorphic in the lower 
complex half plane. 

To state the pole condition for the two dimensional case we map the exterior domain f2 ext C R 2 onto fl v ^ as 
depicted in Figure 1. Here we assume that the material properties are constant on each segment Qj but may vary 
from segment to segment. The z - coordinate remains unchanged under the transformation. The transformed 
Maxwell's equations are exactly of the form (1) but with transformed tensors e^^ and ^ {. With the usual 
notation for the pulled back differential form the weak waveguide problem with transformed exterior domain 
now reads 

Problem 2 (Weak Waveguide Problem with Transformed Exterior Domain). Find k z such that 
there exist fields e(x,y,z) £ i?Q int (curl) and e*(rj,£, z) £ H\oc,k z (curl) such that: 

1. (e*)* = e on the boundary dfl. (Matching Condition) 

2. e*(r?, s) = £e*(?7, £) defines a holomorphic function on the lower complex half plane (3s <0). (Pole Con- 
dition) 

3. The field composed of e and e„* satisfies Maxwell's equations: 

/ (/J,~ 1 die A div - w 2 (ee) A v) + / ( fi~\d\e* A div* - w 2 (e^ ^e*) A v* ) = 

J (ee)Ad Q p+ (e^e.AdopJ = 

J flint J^n.S 



for any v G Hn ititt k z (curl), p £ #n int ,fc a ; an d compactly supported v* £ H\ OCt k z (cml), p* £ k suc h that 
(v*)* = v, (p*)* = p on the boundary £ = 0. 



4. TRANSPARENT BOUNDARY CONDITIONS 



Problem 2 is still posed on an unbounded domain and therefore numerically not feasible. As mentioned in the 
introduction to this paper the transformed exterior field is typically not decreasing in the exterior domain. 
This rules out a simple truncation of the computational domain. When constructing transparent boundary con- 
ditions the aim is to compute the true solution in the interior domain with a numerical effort proportional to the 
number of unknowns in the interior domain. As shown by Schmidt et al. 4 ' 8 the Laplace transform behaves 
very kindly. As numerically approved, a discretization of e* along the real axis with global functions gives a 
transparent boundary condition so that the computed interior solution converges exponentially fast to the true 
solution (up to the interior discretization error) with the number of discretization "points" used for e* . 



In this paper we focus on the Perfectly Matched Layer method introduced by Bcrcnger. 9-11 To motivate the 
method we go back to the one dimensional Helmholtz equation (4). The general solution in the exterior domain 
I + is holomorphic in x. We see that along the straight line (1 + ia) the outgoing part becomes exponentially 
decreasing as far as a is chosen such that a\Rey 'uj 2 pe — k^\ > \^s^uj 2 pe — k 2 \ while the incoming field explodes, 



j£ gi-y/a) 2 /je — k 2 (l+ia)x -L]3 g~ *V w2 M e — k 2 (l+ia)x 

V _ V v ' V v ' 

outgoing ~ evanescent incoming ~ exploding 

Imposing now a zero Dirichlet boundary condition at x — p and assuming that the field intensity is equal to one 
for x = yields 



\B\ 



a i^J u> 2 pe — k 2 (l+i<j)x 



-3?-^/ u 2 [ie — k 2 ap 



g*\/ ui 2 fj,e — k 2 (l+i<r)x _|_ g — *\/ w 2 £ie — k 2 (l+ia)x 

Therefore, the true boundary condition B = is enforced exponentially fast with the layer thickness p. 

In order to switch to the higher dimensional case we assume that e*(?7,<!;) possesses a holomorphic extension 
in £. For a homogeneous exterior domain and some special inhomogeneous exterior domains this is proved in 
Hohage et al.. 7 It is an aim for the future work of the authors to prove that in general a field e„ satisfying 
the pole condition also has an holomorphic extension in £. For 7 = 1 + ia let us denote e*^^, £) = 6,(77,7^), 
£?j,£,b(?7, £) = trj^iriilC), and fi v ^,B(v^ = ^v^iVtlO- The holomorphic extension e^B^, £) is called Berenger 
function. One expects that the field e^B^, £) decays exponentially fast for £ — ► 00. Again this admits to truncate 
the computational domain to f2pML = [?7min, ?7max] x [0,p) and to impose a zero Dirichlet boundary condition at 
£ = p. We are lead to the following PML problem where dk b denotes the exterior derivative with replaced 
by 1/(1 + ia)d s 

Problem 3 (Weak Waveguide Problem with PML). Find k z such that there exist fields e(x,y,z) G 
^fo int ,fcz(curi) and e^PMiX^, £, z) G i?npML,fe z (curl) such that (e*)* — b on the boundary <9f2 (Matching Condi- 
tion) and 

I (p~ 1 d 1 e Adiv- uj 2 (ee) A v) +7 / (\prfi,Be* ; PML A rfi?, - w 2 (e,^3e, jPM L) A v*)) = 

I (ee)AdoP + 7/ (e^Be^pML A d ,BP*) = 

e *,PMLlc=p = 

for any v G Hn iM ,k z (curl), p G H^ iBt kz and v* G i?n PMLife ^(curl), p* G #n PMLife such that (v.)* = v, (p*)* = P 
on the boundary dQ- m t. 

Remark 1. The complex continuation along the straight line 7^ yields a jump in the Neumann boundary 
condition at £ = 0, 

/ (i„| B rf iB e *.BAv»=7 / p,~\d\e* A v*. 




Figure 3. Relative error ||u — lift. 1 12/ 1 1 2 versus thickness of the PML-layer for linear and quadratic finite elements in 
the first experiment. 



The factor 7 left of the integral symbols ^ in Problem 3 is introduced to incorporate this jump in the variational 
problem as the natural boundary condition on dfl m t ■ This avoids the definition of further unknowns on the 
boundary (Lagrange parameters). 

The PML method is proved to converge exponentially fast to the true solution with an increasing layer 
thickness p for a homogeneous exterior domain 12, 13 and for some special inhomogeneous exterior domains. 7 
To demonstrate the accuracy and the exponential convergence of the method even for rather complex exterior 
domains we want to compute the propagation of a TM polarized fundamental mode, 

-AE z -k 2 {x,y)E z =0 

along a waveguide as depicted in Figure 2, see also Zschiedrich et al. 14 The fundamental mode is used as an 
incoming field and is only specified along the left and upper side of the computational domain. Thus this 
example is a non trivial scattering problem - we must recover the propagating mode in the interior domain. The 
exterior domain is non-homogeneous due to the infinite waveguide. Figure 3 shows the relative L 2 - error in the 
computational domain. We observe exponential convergence for growing thickness p of the PML layer until the 
discretization error of the interior problem dominates the overall error. 



5. FINITE ELEMENT DISCRETIZATION 

To discretize Problem 3 we split the interior field 

e = e x (x, y)e ik * z dx + e y (x, y)e lk * z dy + e z (x, y)e lk * z dz 

into a transversal part e_L = e x (x, y)dx+e y (x, y)dy and a longitudinal part e z — e z (x, y)dz. As usual we discretize 
ej_ with Nedelec's edge elements and e z with standard scalar elements. This gives a discrete counterpart to the 
de Rham complex and hence leads to a discrete divergence condition. 15 In this way, spurious modes which may 
rise from the kernel of the Vx - operator when using an improper discretization scheme are ruled out. The 
variational problem for the interior problem reads in classical notation 

f _i [ \7E Z ~ ik z V ± 
J R2 " [ VxxEx 

for all vj_ e Ha iDt (curl) and v z <G #n int with operators 

V± : < nt - #n int (curl) 

(j> (d x cj),dy(j)), 



Vf * — ik z v* 



Lu 2 e 



J 



' E ± " 




' Vv* z 


E z 




ik z v* 



dxdy 
dxdy 



0, 

0, 



Air 




Gold layer 
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Figure 4. Plasmon-Polariton- Waveguide. In the computations we have used a = lOnm, w = 20/im, d\ — 4/im and 
d,2 — 8.01/im. 



V ± x : tfn int (curl) - #Q int (div) 

Ej_ i ^ (d x Ey — d y E x ), 

and 

V±- : #n int (div) - L 2 

D_l i ^ (d x D x + d y D y ), 

Again one sees that any solution to the first equation also solves the second one (divergence condition). Simply 
set vj_ = l/(ikz)Vv z for k z ^ and recall that Vi x Vx = 0. For k z — set v z = and v± = X7p for any 
p G -ffn int ' Within the PML layer we use corresponding finite elements on quadrilaterals, which are defined on 
a reference quadrilateral via a tensor product ansatz. 14 On the whole transformed exterior domain Qpml we 
use a fixed discretization in £ - direction. For the interior discretization we have implemented an adaptive grid 
refinement steered by a residual based error estimator as in Heuveline and Rannachcr. 16 



6. EXAMPLES 

We now demonstrate the ability of our code to cope with challenging problems in the optical waveguide design. 
In the examples, the quadratric eigenvalue problem is solved with the ARPACK package by Sorensen et al. 17 
after a reduction to a linear eigenvalue problem. The Arnoldi method is used in the shift-invert mode and we 
rely on Intel's Math Kernel Library for sparse LU decomposition (PARDISO 18 ). 

6.1. Plasmon Polariton Mode 

As shown by Bcrini ct al. 19 and Bozhevolnyi 20 a very thin metal stripe may serve as a waveguide. In this case 
the propagating mode is localized near the metal stripe. The present geometry is sketched in Figure 4. Since 
the substrate has a relatively high refractive index the modes are typically leaky. Further there are singularities 
near the metal's corner. This calls for an adaptive grid refinement. The coarse grid consists of 13684 triangles 
and is adaptively refined three times during the program execution. Within the PML layer we have used the 
discretization £ = [0.0 : 0.1 : 2.0]. "3 (in Matlab notation). As the initial guess for k z we have used the result from 
the one dimensional problem which is given by a cut along the symmetry axis of the waveguide. In Table 1 the 
computed effective refractive index for the fundamental mode and the computation effort are given. We observe 
convergence up to eight digits after three grid refinement steps. In Figures 5 and 6 one sees isoline-plots for the 
magnetic field strength which show that there are no spurious reflections from the boundary of the computational 
domain. 
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Figure 5. Plasmon-Polariton- Waveguide. Real and imaginary parts of the H z - component for the fundamental mode. 
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Figure 6. Plasmon-Polariton- Waveguide. Left: Absolute value of H z - component for the fundamental mode. There 
appear no spurious reflections at the lower boundary. Right: Magnetic field intensity. The mode is localized near the 
metal stripe. 



Step n PM L,cS N° DOF total time [min] Memory [GByte] 



1.5350262e+00+0.0000981e+00i 159729 01:57 - 0.9 

1 1.5350261e+00+0.0000985e+00i 281151 03:35 - 2.6 

2 1.5350263e+00+0.0000984e+00i 527656 07:52 - 4.8 

3 1.5350263e+00+0.0000984e+00i 881016 12:16 -9.1 



Table 1. Fundamental mode of the Plasmon Polariton waveguide for a vacuum wavelength of Ao = 1.55/im. The 
computations were performed on an AMD Opteron Linux-PC. 



Figure 7. Fundamental leaky mode of the studied Arrow waveguide: Magnitude of the electric field, \E(x,y)\ in the 
cross section. The gray colormap in the right part of the figure is lighted up so that the field in the Arrow layers is better 
visible. Recall that the normal component of the electric field jumps across material boundaries. 



Step 


n>PML,e« 


N° DOF 


n>PML,efi 


N° DOF 





9.9325021c-01+0.0012272c-01i 


51111 


9.9325021e-01+0.0012272e-01i 


51111 


1 


9.9322697e-01+0.0017419e-01i 


92260 


9.9322724e-01+0.0017697e-01i 


135625 


2 


9.9320708e-01+0.0016724e-01i 


154747 


9.9320699e-01+0.0017118e-01i 


404865 


3 


9.9320222e-01+0.0016547e-01i 


265375 


9.9320499e-01+0.0016816e-01i 


1344193 


4 


9.9320574e-01+0.0016710e-01i 


478785 






5 


9.9320580e-01+0.0016820e-01i 


1449444 







Table 2. Fundamental leaky mode of the Arrow waveguide for a vacuum wavelength of Ao = 785nm. The left part 
corresponds to an adaptive grid refinement, the right part to a uniform grid refinement. The computations were performed 
on an AMD Opteron Linux-PC. Computation time and memory requirements are similar to the previous example for a 
equal number of unknowns. Observe that with an adaptive refinement strategy the memory requirements remain below 
the 32-bit PC limit up to the third adaptive refinement step. 
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Figure 8. Hollow core ARROW waveguide. The core width is equal to 12/jm and the core height is equal to 3.5/ira. 
The Arrow layers are composed of silicon nitride (ru = 2.1, w\ = 109nm) and silicon oxide (712 = 1.46, 11)2 = 184nm) the 
substrate has a refractive index of n — 3.4975. 



6.2. Arrow Waveguide 

The present waveguide structure consists of a hollow, rectangular core investigated in Yin et al.. 21 The field is 
confined by antiresonant, reflecting optical layers (ARROW). The geometry is sketched in Figure 8. Again as 
an initial guess we have used the results from the corresponding one dimensional problem on the cut along the 
symmetry axis of the waveguide. Interestingly without transparent boundary conditions we were not able to find 
the two dimensional fundamental mode with primarily TE-polarization. Figure 7 shows the magnitude of the 
fundamental mode. In Table 2 the computed effective refractive index for the fundamental mode is given. The 
adaptive grid refinement allows to compute the propagation mode with a reasonable accuracy even on a 32-bit 
PC. 
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